Time-dependent tunneling of Bose-Einstein condensates 
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The influence of atomic interactions on time-dependent tunneling processes of Bose-Einstein con- 
densates is investigated. In a variety of contexts the relevant condensate dynamics can be described 
by a Landau-Zener equation modified by the appearance of nonlinear contributions. Based on this 
equation it is discussed how the interactions modify the tunneling probability. In particular, it is 
shown that for certain parameter values, due to a nonlinear hysteresis effect, complete adiabatic 
population transfer is impossible however slowly the resonance is crossed. The results also indicate 
that the interactions can cause significant increase as well as decrease of tunneling probabilities which 
should be observable in currently feasible experiments. 
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I. INTRODUCTION 

Since the first experimental preparation of Bose- 
Einstein condensates (BECs) in dilute atomic gases the 
study of their dynamical properties has been a very ac- 
tive field of research. It has become apparent that atomic 
interactions play a crucial role in the prediction or ex- 
planation of a wide range of observable phenomena in- 
cluding, e.g., free condensate expansion, collective excita- 
tions, nonlinear atom optics, solitons, and vortices 
The recent experiment of Ref. || has drawn attention 
to a further dynamical process, namely time-dependent 
tunneling. This work investigated the dynamics of a BEC 
which was accelerated by gravity in the periodic potential 
formed by two counterpropagating vertical laser beams. 
In this way, Bloch oscillations of the condensate were in- 
duced, and each time a turning point of the oscillation 
was reached a fraction of the atoms tunneled into a con- 
tinuum Bloch band. The regular output of atom pulses 
spectacularly proved the macroscopic coherence of the 
initially prepared condensate. 

In the experiment of Ref. pj the influence of atomic 
interactions mainly showed up as a degradation in the 
interference when condensates with high densities were 
studied. As a further consequence, long-time dephas- 
ing of the pulse output was predicted. However, another 
interesting question arising in this context concerns the 
problem of how the atomic interactions affect the indi- 
vidual tunneling processes which arc fundamental to the 
dynamics of the system. The purpose of this paper is 
to work out essential aspects of this question by study- 
ing the modification of the tunneling probability in the 
single process. To this end, the condensate dynamics is 
modeled in terms of a mean-field description provided by 
the Gross-Pitaevskii equation. If the periodic potential 
is sufficiently weak and its period short compared to the 
condensate extension this equation can be simplified by 
means of a two-mode expansion. Thereby, the two modes 
represent the original and the Bragg-scattered contribu- 
tion. In this way one arrives at a set of equations similar 
to the familiar Landau-Zener problem Q. In the present 
case, however, the equations also contain nonlinear terms 



that characterize the effect of the interactions. It should 
be emphasized that these equations, which are central to 
our investigation, are applicable not only to the descrip- 
tion of Bloch band tunneling but are of broader scope. 
They can also model processes such as population trans- 
fer between different hyperfine states by variable exter- 
nal fields ||] or the motion of BECs on coupled potential 
surfaces (which might be of interest in practical applica- 
tions). These processes, as well as Bloch band tunneling, 
are examples of coherent output coupling from BECs (q| 
so that the present work also relates to this context. Fur- 
thermore, it extends recent studies of nonlinear Joseph- 
son oscillations that are based on the time-independent 
version of these equations |i| . The time- independent 
equations are relevant not only in the context of Bosc 
condensation but apply to other areas of physics as well, 
e.g., polaron dynamics, nonlinear optics, biophysics, and 
molecular physics p,|lQ[ so that the present study might 
also be of interest in some of these fields. 

The paper is organized as follows. In Sec. II the model 
is set up and the 'nonlinear Landau-Zener equations' 
mentioned above are derived. Section III first gives a 
brief qualitative discussion of how the nonlinearity af- 
fects the tunneling probability. We then discuss the main 
result of our investigation which shows that for certain 
parameter values the tunneling probability does not van- 
ish however slowly the system evolves. In other words, 
a complete adiabatic population transfer is impossible 
under such circumstances. This effect which would not 
be possible in linear two-level crossing models is a direct 
consequence of the discrete self-trapping transition [jlO| . 
Using a phase space representation of the problem the 
tunneling probability for very slow processes can be de- 
termined. Subsequently, we discuss the population trans- 
fer for fast and slow resonance crossing using numerical 
methods and simple analytical models. In Sec. IV the 
predictions of the two-mode system are compared to the 
numerical solution of the Gross-Pitaevskii equation de- 
scribing Bloch-band tunneling. These studies also show 
that significant modifications of the tunneling probability 
due to nonlinear effects should be observable in currently 
feasible experiments. Conclusions are given in Sec. V. 
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II. THE MODEL 

The Gross-Pitaevskii equation for a condensate wave 
function tp(r,t) undergoing Bloch band tunneling in a 
periodic optical potential is given by 



ihip = Hi in 4> + g\ip\ 2 ip 



with (lljjl 



Hi, 



h 2 v 2 

2m 



+ (V T + V opt -Fz). 



(1) 



(2) 



The atomic mass is denoted m, and g — Ai:Ti 2 aN/m is 
the nonlinearity parameter with a the s-wave scatter- 
ing length and N the number of atoms in the sample. 
The condensate is assumed to be tightly confined in the 
radial direction by the trapping potential Vt, the z de- 
pendence of which is neglected. The periodic optical po- 
tential is given by V opt = Vq cos(2fc[,z), with the laser 
wave vector Ul , whereas the potential inducing the Bloch 
oscillations is characterized by the accelerating force F. 
Such a potential may be produced by tilting the stand- 
ing wave || or by frequency-shifting the counterpropa- 
gating laser fields jljj. We want to discuss the conden- 
sate time evolution under the following two conditions. 
First, the condensate has an initial momentum pq well 
defined on the scale set by the lattice wave vector 2fc_L, 
i.e., l z 3> A = n/ltL with l z the axial extension of the 
BEC. The time evolution starts with the BEC well sepa- 
rated from any tunneling resonances, i.e., avoided cross- 
ings between Bloch bands. These resonances occur for 
condensate momenta around (21 + l)Uki,, I = 0, ±1, . . . 
For concreteness we take po = in the following [apart 
from the example of Fig. ||(b)] ■ The condensate is studied 
as it passes through the resonance at Kkh- The second 
condition is that the optical potential is sufficiently weak, 
i.e., Vq <C Ti 2 k 2 L /2m. In this case, the crossing of the reso- 
nance can be described in terms of a coupling to a single 
state with a momentum shifted by —2hkL. The wave 
function ip(r,t) is expanded in terms of the original and 
the Bragg-scattered contribution, i.e., 



Vj = e i{kL+Ft/n) Z(j)+ 

+ e i(-kL+Ft/h)z 



hkr 



Ft' 



(z-—t-^—)b + (t) (3) 



( Tik L Ft 2 \ l .. 
(z + —t-—)b_(t) 
V m 2m / 



with t = corresponding to the point of exact reso- 
nance and b+ = 1, 6_ = 0, initially. The two enve- 
lope functions <f>± (whose radial motion is frozen due to 
the tight confinement) are normalized to one. It is as- 
sumed that 0_ is very similar to cj>+ , the shift in position 
between the two being negligible at all times t relevant 
for the tunneling process, i.e., we set J d 3 r<f)-(z)(j)+(z — 
2hkLt/m) ~ 1. This approach is valid if the nonlin- 
ear effects are not too large and the tunneling process 
is sufficiently short. Expression (||) is then inserted 



into Eq. (|j) thereby discarding second derivatives of <p± 
(slowly- varying envelope approximation). By projecting 
the resulting equation onto the instantaneous wave func- 
tions e l ^ kL+Ft ^ h ^ z 4>± and switching to a rotating frame 
where a±(t) = exp[i J u(t')dt']b±(t), u(t) = fik 2 L /2m + 
F 2 t 2 /2mh — et/2 — 2-f, we obtain the 'nonlinear Landau- 
Zener equations' 



ia + = eta+ + 17a_ + 7|a + | 2 a+, 



with 



2Fk, 



= f2a+ + 7|a_ 



2h 



-| jd\\^{v)\\ 



(4) 



(5) 



In the following, f2 > is assumed. As mentioned in the 
Introduction similar equations can be derived in other 
contexts, e.g., for population transfer between different 
hyperfine ground states with variable external fields [g]. 
In Sec. IV it is shown that Eqs. (Q) allow an accurate pre- 
diction of tunneling probabilities in realistic situations. 

Significant insight into the system behavior can be ob- 
tained by noting that Eqs. (^) may be derived from the 
Hamilton function Il7|,|l0| 



H(N+,@) 



AN+ + 2n^N + {l - N+)cos 6 
+ 7 (Nl -N+ + 1/2) (6) 



with iV+ = |a+| 2 and O — arg(a+al) as canonical vari- 
ables and A — et. The dynamics induced by this Hamil- 
ton function for fixed, time-independent A is discussed 
in detail in Rcfs. ||7|,p|Jl0|. In the present context the 
stationary states S = (Ns,®s) 01 H(N + ,Q) are of par- 
ticular interest because they take over the role of the adi- 
abatic eigenstates in the linear problem (7 = 0). From 
dH/dQ = it follows that they always have Qs = 
or 7T. The condition dH/dN + = then shows that 
Ns = n s + 1/2 is obtained as a real- valued solution of 
the equation 



A 



I 2 



A!_ 

4 7 2 



47 



A 2 
I67 2 



= 0. 



(7) 



If I7I/O < 2 there are exactly two stationary states S- 
and S + having 6 = and n, respectively, for all values of 
A. They correspond to the high- and low-energy eigen- 
states of the linear problem. For |7|/£1 > 2 two further 
states appear in the vicinity of A = as discussed in 
more detail in Sec. III. 



III. TUNNELING PROBABILITIES FOR THE 
NONLINEAR LANDAU-ZENER MODEL 
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A. Overview and qualitative discussion 

In the study of the tunneling probability for the non- 
linear Landau-Zener model we are interested in the long- 
time solution of Eqs. (||) provided the system is prepared 
at t — > — oo in the high- or low-energy (stationary) state. 
We choose the low-energy state a+ in the following, i.e., 
a + (t — > — oo) = 1, and the tunneling probability is thus 
Pt = \a+(t — > oo)| 2 . In the linear case (7 = 0) the 
tunneling probability is independent of whether the sys- 
tem starts in the low- or the high-energy state. In the 
nonlinear case the same still holds true if the sign of the 
nonlinearity is changed as well, i.e., the tunneling prob- 
ability for a system starting in '+' is the same as for one 
initially in '— ' and having a nonlinearity parameter —7. 
The whole range of dynamical phenomena predicted by 
Eqs. §) (for both signs of 7) can thus be observed with 
a condensate of repulsive interactions by using both ini- 
tial conditions. Due to the scaling properties of Eqs. (Q) 
Pt can be considered a function of fi/e 2 and 7/O, only. 
An overview of the dependence of Pt on these parame- 
ters is given in Fig. |l| which shows results obtained from 
the numerical solution of Eqs. (^). In the linear problem 
the numerical calculation of Pt is facilitated by trans- 
forming from a + , a_ onto the basis of the instantaneous 
'dressed' eigenstates. In this basis the oscillations of the 
system trajectory around its asymptotic limit for t — > 00 
are suppressed and after a relatively short propagation 
time the value of Pt can be read off from the population 
of the eigenstates. This approach is generalized to the 
nonlinear case by computing the canonical action of the 
trajectories as determined by the Hamilton function 
This action plays the role of the adiabatic invariant and, 
in the linear case, is equivalent to the eigenstate popula- 
tion. Various aspects of the behavior of Pt are discussed 
below. 

For the linear problem the tunneling probability is 
given exactly by the well-known Landau-Zener formula 
Pt = cxp(— 2-itVL 2 /e) [Q. A qualitative understanding of 
how the nonlinear terms influence Pt can be obtained 
by noting that they give rise to an effective detuning 
A e // = et + 7(2|a + | 2 — 1) between the states a+ and 
a_. Unless the resonance is traversed too rapidly the 
system remains in the vicinity of S+ (st) and the popula- 
tion |a+(i)| 2 closely follows its steady state value Ns + = 
n Sy+ (et) + 1/2. A rough estimate of the tunneling prob- 
ability can be obtained from the rate R — A e ff(t = 0) 
at which the point of zero detuning is crossed. As a first 
approximation it follows from the linear Landau-Zener 
formula that Pt ~ exp(— 2ttQ 2 /R). As for small A 

n St+ m -A/[2(2fi + 7)] + QA 3 /[2(2n + 7 ) 4 ] (8) 

one thus finds 

ft 2 

P T «exp{-27r — (1 + 7/20)}. (9) 




1,2 3 



FIG. 1. Tunneling probability Pt as a function of fi 2 /e for 
various values of 7/fl. 

Although not quantitatively accurate this formula indi- 
cates two trends confirmed by the detailed investigation: 
(i) nonlinear effects become significant as soon as 7 be- 
comes comparable to f2 (as one might expect) and (ii) 
Pt is increased (decreased) for 7 > (7 < 0) as com- 
pared to the linear case. Interestingly, the expansion of 
n Sj + also shows that the third-order correction to A e // 
renders the tunneling transition superlincar (sublinear) 
for 7 > (7 < 0) in the terminology of Ref. ||l4| , i.e., 
the resonance is effectively crossed faster (slower) than 
the linear approximation predicts jl5). The results of 
Ref. |l4| indicate that the approximation (^|) should un- 
derestimate the nonlinear effects on Pt- This is indeed 
observed in the present case if VL/e 2 is not too small and 
7/fi sufficiently larger than —2. 

B. Nonlinear hysteresis effects in the tunneling 
probability 

The above discussion already indicates the significance 
of the stationary state behavior for an understanding of 
the tunneling probability. The breakdown of the expan- 
sion for n Si+ at 7/fi = —2 now suggests that at this point 
a qualitative change in the system behavior occurs. In- 
deed, for j/Q < —2, as A is increased from large negative 
values, two further stationary states of H(N + ,Q), S' + 
and U+, emerge at a detuning — A c < [see Fig. ||(a)]. 
They both have the same phase = n as S+ and pop- 
ulations Ns>,+ < Nu,+ < The point S' + is stable 
whereas U+ is unstable. With growing A, U+ approaches 
S+. At A c they coalesce and only S'< remains (besides 
S-). We thus observe a typical nonlinear hysteresis phe- 
nomenon. How does this scenario affect the tunneling 
process? As long as 7/fi > — 2, for e — > the system will 
stay arbitrarily close to S+ over the whole evolution of 
the system so that Pt — > 0. For 7/0 < —2, however, the 
hysteresis effect gives rise to a different behavior which is 
sketched in Fig. |2|(b). For small e the system will again 
remain close to 5+, initially. However, as U+ approaches 
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FIG. 2. (a) Coordinate N+ of stationary states as a fun- 
cion of A for 7/fi < —2. (b) Schematic phase space plot in 
the vicinity of A c . The system trajectory eventually has to 



switch onto the orbit C. 
a function of — 7/Sl. 



(c) Asymptotic probability P^ 00 ' as 



S+ it eventually has to switch from its trajectory in the 
vicinity of S+ onto a large orbit C encircling S' + [ p^| . 
Subsequently, the system evolution is quasi-adiabatic as 
t — > 00. The final tunneling probability is thus deter- 
mined by the canonical action of the orbit C which is an 
adiabatic invariant in this case. We thus conclude that 
in the limit fi 2 /e — ► 00 the tunneling probability is given 
by 



= P T (&/e -> 00) = J- 

Z7T 



SdN 4 



(10) 



with C* the orbit passing through U+ = 5+ at A = 
A c and encircling S' + . For small e > the orbit C, 
onto which the system switches, completely encloses C* . 
Therefore, P^ 00 '' can be expected to be a (at least local) 

minimum of the tunneling probability. In Fig. ||(c) P T °°^ 
is shown as a function of — 7/fi. It becomes apparent 
that, even for modest values of —7/fi, P T °°^ is not small 
compared to one. The above discussion is confirmed by 
the numerical simulation of Eqs. (^|). In particular, it is 
found that for fixed —7/fi < —2, Pt is a monotonically 
decreasing function of fi 2 /e (cf. Fig. Jib. As fi 2 /e — > 00, 
it tends to the value given by Eq. (10) which is thus a 
global lower bound on the tunneling probability. Further- 
more, for small e, the canonical action along the system 
trajectory abruptly changes around A c , as expected. 



C. Rapid passage through resonance 

Having discussed the asymptotic limit of the tunnel- 
ing probability we now study more closely its behav- 
ior for small values of fi 2 /e, i.e., when the resonance 
is passed rapidly. In this case some insight may be 
gained from a perturbative analysis of Eqs. @) with fi 
the small parameter. In a rotating frame with a + = 



a + exp(— iet 2 /2 — i'ft), a_ = a_ Eqs. (|J) read 
ia + = fia_ exp(ie< 2 /2 + i'jt) + 7(|a+| 2 — l)a + , 



1X1- 



fia+ exp(— iet 2 /2 — i'yt) + 7|a_| 2 a_. 



(11) 
(12) 



These equations are now iterated in the standard way 
with initial conditions a^°'(— 00) = 1, a_ (— 00) = up 
to third order. This yields 



d { ^(t) - -iQ / <fti£(ti) 

J — OO 

rt nil rte 



(13) 



-in- 



—00 j —00 

t rti 



dti / dt 2 / A 3 £(ii)£*(t 2 )£(i 3 ) 



7fi 3 / dh 



dt 2 £(t 2 



fiTi + fi 3 T 3 + 0(fi 5 ) 



2 fti 



dt 2 £(t 2 ) + o(n 5 ) 



with £(t) = cxp(— iet 1 1 2 — ijt). The first integral con- 
verges for t — > 00 whereas the other two diverge. This be- 
havior is to be expected because, e.g., in the linear prob- 
lem (7 = 0), although the moduli of a± converge, their 
phases contain logarithmically divergent terms which are 
reflected in the divergence of the fi 3 contributions. A 
convergent approximation of Pr up to order fi 4 can nev- 

ertheless be obtained from a_ (t — > 00). Writing Pt ~ 

l-|a_ 3 (*^oo)| 2 = l-fi 2 |T 1 | 2 -2fi 4 Re(T*T 3 ) + 0(fi 6 ) 
it is found that the expression Yle(T*T^) converges; all di- 
vergent contributions in the fi 3 integrals are contained in 
\m.(T*Ts). Determining the limits t — > 00 of the relevant 
expressions the approach finally yields 



P T = 1 - 2tt- 



fi 2 



2tt' 



,fi 



,7fi 4 



£2 -^ E , 5 



+ 0(fi 6 ). (14) 



The first three terms are familiar from the expan- 
sion of the Landau-Zener formula Pt — exp(— 27rfi 2 /e) 
for the linear problem. The fourth term repre- 
sents the first nonvanishing contribution of the non- 
linearity. The numerical coefficient J\f is equal 
to Re[e-"/ 4 yi28^ J^dtl^it)] 2 ^^)} « 15.751 with 

At this point it has to be emphasized that it is not 
claimed that Eq. ([l4|) is rigorously valid as expansion 
of Pt around fi = 7 = 0. Nevertheless, it yields the 
following useful information: (i) For very rapid passage 
through resonance the nonlinear interactions do not influ- 
ence the tunneling probability. Their contributions arise 
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FIG. 3. Numerical calculation (full curves) of Pt and 
approximation according to Eq. (h4|) (dashed curves) for 
7/ft = -10.0, 0.0, 5.0. The result of Eq. @ for j/Sl = 10.0 is 
given by the long-dashed curve that is close to the numerical 
one for ~y/Q = 5.0. 



only in higher order in ft (see Fig. H). (ii) Equation (E 
correctly predicts that for 7 > (7 < 0) Pt is dimin- 
shed (increased), (iii) Quantitatively, Eq. (|lj) gives a 
good approximation of Pt in particular for negative 7 
and could be used, e.g., at 7/ft = —10 to estimate Pt 
for ft 2 /e < 0.04 (cf. Fig. ^). For positive 7 the agreement 
is not as good. In the interesting regime 1 < 7/ft < 10, 
ft 2 /e < 0.04, the nonlinear effects are underestimated 
by a factor of about 2. This is because the influence of 
higher-order terms in 7, which are neglected in Eq. (14), 
is more relevant for positive 7. 



D. Slow passage through resonance 

To examine the tunneling probability for slow passage 
through resonance in more detail we again use the idea 
that under these circumstances the system stays close to 
the instantaneous stationary state S+ . A simple analyti- 
cally tractable model is obtained if in Eqs. (0) |a+| 2 is re- 
placed by the stationary state population Ns,+(t) of the 
linear system, i.e., \a + \ 2 « 1/2 - A/[4(A 2 /4 + ft 2 ) 1 / 2 ] 
with A = et. The tunneling probability PyP for this 
system can be estimated with the help of the semiclas- 
sical theory of nonadiabatic transitions flJIO] "H which 
is based on the study of the quasienergy function Q(t) = 
[A 2 eff (t)/4 + Sl 2 ]V 2 . Thereby, A e// = + 7 (2|a + | 2 - 1). 
If tk are the zeroes (or complex crossing points) of Q{t) 
in the upper complex half plane the tunneling probabil- 
ity can be approximated as P T 1 ^ — | exp(i2?fc)| 2 
with T>k = 2 f* k Q(t)dt and = — 2i lim t ^ tfc (t — 
tk)^iA e ff(t)/Q 2 (t). As briefly shown in the Appendix, 



for |7|/0<1, Prp is well approximated by 

= exp(-/')[2cos(27r7ft/3e) - l] 2 , 7 > 0, 
P® = exp(-Z'), 7 < 0, 




FIG. 4. Tunneling probability for j/Sl = ±0.25. Bold 
curves: exact results obtained from Eqs. (^), short-dashed 
curves: approximation according to Eqs. (jlil). The 
long-dashed curves show the tunneling probability with |a+| 2 
in Eqs. (^) replaced by the exact stationary state population 
as determined from Eq. (0). For larger |7|/fi the discrepancy 
is even more pronounced. The thick line shows Pt for 7 = 0. 



with /' = 2f2 2 {7r - 7 [| + 21og(^ r 7|7ft/2)]}/e. Com- 
parison to the solution of Eqs. (Q) shows that the formula 
for 7 > gives a useful approximation to Pt well before 
the first minimum of P^, i.e., if ft 2 /e<ft/47 (see Fig. 
^). For negative 7 Eqs. (|l|) typically underestimate the 
nonlinear effects but they yield a rough first estimate for 
7/Q less than about 0.25. More interestingly, however, 
Eqs. ( |l5| ) show a general behavior which is observed for 
the tunneling probability of the nonlinear problem as well 
(see Fig. |). If 7 < 0, P T l) is a monotonically decreasing 
function of ft 2 je whereas for 7 > it is oscillating. This 
behavior is a consequence of the structure of the com- 
plex crossing points of the quasienergy. For 7 > three 
crossing points contribute to P T thus causing interfer- 
ence effects whereas for 7 < there is only one. With 
a theory similar to the one of Ref. |17| not available for 
nonlinear systems the study of this model problem gives 
an idea of how these features in the behavior of Pt may 
come about and shows that they appear in a broader 
class of systems. Numerical studies of the original Eqs. 
([|) indicate, however, that at fixed 7/ ft > the min- 
ima of Pt approximately occur at positions kCfl 2 je with 
k = 1,2,..., and C a function of 7_/ft- This is in clear 
contrast to the behavior of Eqs. (fi~q), furthermore, the 
first minimum appears at a much larger distance from 
ft 2 /e = than indicated by these equations. One might 
assume that it is possible to improve on the above results 
by studying Eqs. (||) with |a+| 2 replaced by the exact sta- 
tionary state population Ng + (t) as determined from Eq. 
(Q). Unfortunately, the analytical theory of Ref. [jlTj is 
not applicable in this case as no complex crossing points 
appear. The numerically obtained tunneling probabil- 
ity has the same qualitative features as discussed above; 
quantitatively, however, the approximation to the exact 
Pt is not as good as one might intuitively expect (cf. Fig. 
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IV. BLOCH BAND TUNNELING WITH 
BOSE-EINSTEIN CONDENSATES 

In this section the prospects of experimentally observ- 
ing the nonlinear effects discussed in Sec. Ill are exam- 
ined with a focus on Bloch band tunneling. To this end a 
one-dimensional simulation of Eq. ([!]) for a sodium con- 
densate is studied. For the calculation, a condensate ex- 
tension of l z — 160/im and a period of the optical poten- 
tial A = 0.5/im are used. With £1 = 10 4 s _1 higher-order 
momentum components are detuned by an amount of at 
least 40O from the '+/— ' modes when the resonance is 
crossed. The two conditions for the applicability of ap- 
proximations (||), (^J) which are given in Sec. II are thus 
well fulfilled. When the wave packet reaches the tunnel- 
ing resonance at 3hkT that follows the initial one at Tik L 
the two momentum components + and — have a detun- 
ing of 40f2, i.e., the tunneling process is well defined. For 
il 2 je w 1 the acceleration is of the order 10ms -2 . The nu- 
merical simulations show that on the time scale of passing 
through the resonance (less than 10~ 2 s in the examples) 
the spreading of the initial wave packet as well as the 
spatial shift between the two modes is small. The crucial 
parameter 7 takes on the value 10 _16 n c m 3 s _1 with n c the 
condensate density. In our example a ratio of 7/fi = 1 
can be reached for a BEC with N — 5 x 10 6 atoms and a 
radius of 10/im. All the parameter values given are well 
within the realm of currently feasible experiments . 

The numerical simulations of Eq. (]]]) are performed 
with a standard split-step algorithm (see, e.g., [pOfl ). 
Thereby, the initial state is taken as the ground state 
of a harmonic potential yielding the required value of l z . 
The main results of these calculations are shown in Figs. 
||. In these diagrams the fractional population P + of the 
initial momentum component is displayed as a function 
of r = tme/Ah 2 k L1 i.e., time measured in units of the 
temporal distance between two resonance crossings. The 
point t = corresponds to zero detuning between the 
two relevant momentum modes. Numerically, P+{t) is 
determined from the Fourier transform of tp(z,t). Ex- 
perimentally, it could be obtained by switching the pe- 
riodic potential off at t and observing the subsequent 
condensate evolution. In Fig. @(a) condensates with var- 
ious values of 7/fi are started at rest with the acceler- 
ation chosen such that the standard Landau-Zener for- 
mula (disregarding the atomic interactions) predicts a 
tunneling probability of 0.01. Having passed the reso- 
nance -P+ quickly settles down to oscillations around a 
well-defined mean value which may be regarded as the 
tunneling probability. These mean values depend some- 
what sensitively on the initial conditions. For exam- 
ple, they may change by up to 0.02 if the initial veloc- 
ity is varied by zL0.01hkL/m. Nevertheless, they com- 
pare reasonably well to the asymptotic tunneling prob- 
ability Pt (indicated in Fig. || with filled circles •) ob- 
tained from Eqs. (^) for the corresponding parameter val- 
ues. Typically, the values for Pt are smaller than the 
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FIG. 5. Fractional population P+ of the initial momen- 
tum component as a function of scaled time r = tme / '4h 2 'k\ 
for various values of 7/fl. In (a) the initial wave packet 
velocity u, = and the acceleration is chosen to obtain a 
tunneling probability of Pl = 0.01 according to the Lan- 
dau-Zener formula in the absence of nonlinear interactions. 
In (b) Vi = —2h/kLm and Pl = 0.1. The predictions of Eqs. 
(0) for Pt at the respective parameter values are indicated by 
filled circles •. 

results obtained from the Gross-Pitaevskii equation by 
approximately 0.01—0.02. Most importantly, however, it 
is seen that for 7/fi = 1.5 the tunneling probability is in- 
creased by almost an order of magnitude compared to the 
case 7/fi = 0.01. These results show that Pt is very sen- 
sitive to atomic interactions in experimentally relevant 
situations. Figure ||(b) illustrates that one can also ob- 
serve a significant decrease in the tunneling probability. 
As discussed in Sec. III. A, this occurs if the condensate 
starts from the high-energy state, i.e., in the present sit- 
uation it needs to be launched with a velocity less than 
— ^fci/m. In Fig. ||(b) the initial velocity —2fik L /m was 
used which places the wave packet in the middle of two 
tunneling resonances. The acceleration was chosen such 
that without nonlinear interactions a tunneling probabil- 
ity of 0.1 is expected. For smaller values of 7/O there 
is again good agreement between the predictions of the 
two-mode model and the solution of the Gross-Pitaevskii 
equation. At 7/fi = 2.0, however, Eqs. (0) predict a tun- 
neling probability of about 10 -4 whereas the simulation 
yields a significantly higher value of 0.01. 
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V. CONCLUSION 

In this paper the influence of atomic interactions on 
time-dependent tunneling processes of Bose-Einstein con- 
densates was investigated. A nonlinear Landau-Zcncr 
equation was derived which describes main aspects of 
processes such as Bloch band tunneling, ground state 
population transfer with variable external fields and 
condensate motion on coupled potential surfaces. The 
tunneling probabilities predicted by this model were 
discussed in detail, in particular, it was shown that 
for strong enough nonlinearities a complete population 
transfer by adiabatic following is impossible. This behav- 
ior is a consequence of a nonlinear hysteresis effect. To as- 
sess the reliability of the nonlinear Landau-Zener model a 
comparison to simulations of the Gross-Pitaevskii equa- 
tion was made. 

In actual experiments it is only possible to determine 
the tunneling probability for a certain portion of the 
(fl 2 /e, parameter space. For example, the accelera- 
tion cannot be made arbitrarily small so that the asymp- 
totic behavior of the tunneling probability shown in Fig. 
^(b) may never be verified directly. Furthermore, tun- 
neling probabilities close to or 1 are very difficult to 
measure accurately due to the long-time oscillations of 
the mode populations. Nevertheless, as shown in Sec. 
IV, drastic nonlinear effects already appear in the exper- 
imentally accessible parameter space and they are appro- 
priately described with the simple model of Eqs. (^). Its 
detailed study is thus well justified. 

However, to obtain a full understanding of time- 
dependent tunneling it is necessary to expand the in- 
vestigation beyond the limits of applicability of Eqs. @). 
In the context of Bloch band tunneling new features in 
the system behavior may arise, for example, if the exten- 
sion of the condensate becomes comparable to the period 
of the optical potential. Another potentially interesting 
question concerns the study of quantum mechanical ef- 
fects beyond the mean field description that was applied 
here. 

This work was supported by the United Kingdom En- 
gineering and Physical Sciences Research Council. 



The integral T>i is approximated by expanding the inte- 
grand as 



APPENDIX: DERIVATION OF EQS. (15) 



The complex crossing points tk are obtained as solu- 
tions of 



W 4 \ 2^A 2 /4 + ft 2 



0, (16) 



with A = et, which can be converted into a quartic equa- 
tion. If 7 < there is only one root in the upper complex 
half plane which is given by 

A 1 = ^[l-(| 7 |/0) 2 / 3 /2 + 0(|7| 4/3 )]. 



Q(t) « VA 2 /4 + f> 2 - 7 A 2 /[8(A 2 /4 + n*)]. 

The integration can then be carried out analytically and 
the second of Eqs. ( [l5| ) is obtained with Ti = —1. 

In the case 7 > Eq. ( |l6| ) has three solutions in the 
upper complex half plane. They are given by 

Ai, 2 = 2iil{l + e ±M / 3 ( 7 /fi) 2 / 3 /2 + 0( 7 4 / 3 )}, 



A 3 = 2iVL 



1 - 



if 
32 ft 2 



0( 7 4 ) 



with Ti.2 = —1 and Tj, = 1. The integrals P1.2 are 
determined as in the case 7 < 0. The integral T>% cannot 
be dealt with in this way. However, a numerical analysis 
shows that for 7/Q < 1 one can set to a good degree of 
approximation T>i = ilm(2?i 2)' This approach yields the 
first of Eqs. (fi"5|). The comparison of these equations to 

the numerically determined p!p confirms the accuracy of 
the approximation for |7|/fi < 1. 
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